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We present simulations of binary black holes mergers in which, after the common outer horizon has 
formed, the marginally outer trapped surfaces (MOTSs) corresponding to the individual black holes 
continue to approach and eventually penetrate each other. This has very interesting consequences 
according to recent results in the theory of MOTSs. Uniqueness and stability theorems imply that 
two MOTSs which touch with a common outer normal must be identical. This suggests a possible 
dramatic consequence of the collision between a small and large black hole. If the penetration were 
to continue to completion then the two MOTSs would have to coalesce, by some combination of the 
small one growing and the big one shrinking. Here we explore the relationship between theory and 
numerical simulations, in which a small black hole has halfway penetrated a large one. 

PACS numbers: PACS number(s); 04.20Ex, 04.25Dm, 04.25Nx, 04.70Bw 



2 


I. INTRODUCTION 


An unexpected feature has been observed in the simulation of equal mass binary black holes following their inspiral 
and merger [1]. After formation of a common apparent horizon, as the two marginally outer trapped surfaces (MOTSs) 
corresponding to the individual black holes continued to approach, they eventually touched and then penetrated each 
other. This penetration was surprising because it had not been considered in theoretical discussions and had not been 
observed in prior simulations of binary black hole mergers. In retrospect, it is tempting to speculate in some heuristic 
sense that a small black hole should enter a very large black hole with hardly any notice of its presence. In fact, 
it has been been conjectured on the basis of the equivalence principle that a very small black hole should, in some 
appropriate sense, fall into the large one along a geodesic. However, such a perturbative picture is unreliable in the 
interior of the event horizon surrounding the two black holes, where the MOTSs exist. Here we use an independent 
evolution code to first confirm that the individual MOTSs do penetrate following a binary black hole merger and we 
analyze some of the highly interesting features revealed by the simulations. 

The simulations which first demonstrated the penetration of the individual MOTSs were carried out with a code 
based upon Fock’s treatment [2] of the harmonic formulation in terms of the densitized metric i.e., the PITT 

Abigel code which had been developed for the purpose of studying outer boundary conditions 011]. By incorporating 
adaptive mesh refinement and a horizon tracker available in the CACTUS toolkit |S], the code was further developed 
to simulate a binary black hole inspiral using excision to deal with the internal singularities. The simulations presented 
in the present paper are based upon an independent harmonic code using standard 3 + 1 variables and treating the 
singularities via punctures [5| rather than excision |7]. We describe the computational methods in Sec. [iTj 

The seminal work of Hayward 0 and Ashtekar and Krishnan 0 has led to a rich mathematical theory of the 
dynamical horizons traced out by the evolution of MOTSs, as reviewed in [lOj . This theory has provided insight into 
both the classical and quantum properties of black holes. In particular, the subject has strong bearing on numerical 
relativity because a MOTS can be identified quasilocally by the vanishing expansion of its outward directed null rays, 
whereas the identification of an event horizon requires global information which is not available at the early stages of 
a simulation. Thus MOTSs play a central role in the preparation of initial data sets for black holes and in tracking 
their subsequent evolution. A numerical study of MOTSs has been carried out nmn], which confirmed the main 
features expected from the theory regarding the early stages of a binary inspiral. Recently, there has been substantial 
new theoretical development centered around the uniqueness and stability of MOTSs [T5HT5] . This recent theory 
applies to the general case of the marginally outer trapped tube (MOTT) traced out by a MOTS, which need not 
be spacelike as in the special case of a dynamical horizon. It has important bearing on how MOTSs approach each 
other and penetrate. We review the main mathematical results and their relevance to the binary black hole problem 
in Sec. IHIl 

The simulations presented in Sec. IV confirm the results observed in |T] and significantly extend the degree of 
penetration of the MOTSs. There are five distinct stages as the MOTSs approach and penetrate. 


• The large separation of the individual MOTSs. The properties of this stage are mainly determined by the choice 
of binary black hole initial data. 

• Formation of a common apparent horizon as the MOTSs approach. 

• The moment of external osculation of the two MOTSs as they touch. 

• The penetration of the two MOTSs. 

• The ultimate fate of the individual MOTSs. 


The original simulation showing the penetration of two inspiraling MOTSs [T] and the numerical study of MOTSs 
in the early stages of a binary black hole [T^] were confined to the case of equal mass black holes. Here we consider the 
unequal mass case, with the goal of shedding light on what happens in the extreme mass ratio regime. Because of the 
computational cost in resolving different length scales, we limit ourselves to the case of a mass ratio rngrnaii/i^iarge = 
1/4. In order to better understand the geometry of the merger, we also restrict ourselves to the simplest case of a 
head-on collision. Already in this case, there are quite complicated geometrical effects. 

Our simulations are based upon time symmetric initial data for the two black holes. The time symmetry introduces 
some non-intuitive features of the initial MOTSs, particularly when they are close together, which emphasize the 
importance of looking at their invariant geometrical properties rather than their coordinate description. We repeat, 
for the unequal mass case, the simulation of the head-on collision carried out by Krishnan and Schnetter |12| for equal 
mass black holes. They simulated this equal mass case with a code based upon a 3 + 1 formulation of Einstein’s 
equations (the AEI BSSN formulation), which is quite different from the harmonic formulation used here. As a check 
on our code, we find qualitative agreement with the results of [H] for the early stages before the individual MOTSs 
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touch. We then continue the simulation of the head-on collision to the later stage where the individual MOTSs 
penetrate, which could not be treated in [12]. 

The original simulations of the penetration in [T| were limited by the excision region inside the individual MOTSs. 
The present code is able to track the penetration to a much further stage, essentially until the small MOTS has 
penetrated halfway into the large MOTS. At that stage its interior puncture is close to the surface of the large MOTS 
and the horizon finder is not able to track it. The simulation of the penetration is validated by the convergence 
measurements given in Sec. |IVE[ 

We demonstrate by numerical simulation that before the individual MOTSs make contact a common outer apparent 
horizon forms, in accord with the theory described in Sec. [ml When the two MOTSs first touch, the theory also predicts 
that their mean surface curvatures must agree at the point of osculation. Consistent with this theory, our simulations 
show that the small MOTS creates a strong distortion of the large MOTS at the contact point which causes the mean 
curvature of the large MOTS to grow rapidly at the point of contact, i.e. its curvature radius rapidly shrinks. On 
the other hand, at this stage the large MOTS has only a small effect on the small one. We track this growth of 
mean curvature of the large MOTS as the penetration proceeds. In a scenario where the black holes are produced by 
collapsing matter, without any puncture, one might expect the penetration to continue to the point where the small 
MOTS completely penetrates the large one. If this were to occur, then a dramatic consequence of the underlying 
uniqueness theorems for MOTSs would ensue. At the time where the back of the small MOTS enters the front of the 
large one, the theorems imply that the two MOTSs must identically coalesce, either through shrinkage of the large 
MOTS or growth of the small MOTS. Although our simulations cannot proceed to this stage, they provide some 
glimpse of how it might proceed. 


II. COMPUTATIONAL METHOD 
A. The evolution system 


We employ an evolution system based upon the 3-1-1 formulation of Einstein equations in generalized harmonic 
coordinates = {t, x*) as described by Friedrich and Rendall [T7]. (We make minor alterations in the notation of [T7j 
to conform to standard usage in numerical relativity.) This 3-1-1 foliation gives rise to the metric decomposition 



(2.1) 

where 


= adfj,t, a = — 

V-5 

(2.2) 

is the unit future directed timelike normal to the foliation. The evolution proceeds along the streamlines of the vector 
field dt = anf^dfj^ + determined by the lapse a and shift = (0, /?*). 

In this formulation the constraints are 

Qfi 

(2.3) 

where E^ is related to the Christoffel symbols by T^ = and F^{g, x) are harmonic gauge 

When the harmonic constraints are satisfied, i.e. = 0, the Einstein equations reduce to a 

wave equations for the spatial metric hij, shift /3* and lapse a which take the form |17j 

source functions [18] . 
system of quasilinear 

- g^^'df.dj.hij = Sij 

(2.4) 

^ (dt - 

(2.5) 

^ (dt - P^djf a - DjD^a = S, 

(2.6) 


where Di is the covariant derivative associated with hij and the right-hand-side S-terms do not enter the principal 
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part. In detail, these terms are 

2 2 

Sij = -^Kij {dta - Cpa) + -^DiaDja (2.7) 

-2Z7(, [h,)khtF’'] + {dtP'^ - P^deP'^) 

+ ^ ( 5 ,)/ 3 ^) (d,P^) 

) ^I3hi)k — {^(jP ) {9ehi)k) 

+4K,,K^ - 2K,,K - 


S'* = 4 {K^‘^ - KF^) Dia 

-2a - Kh^^) - {d,P^) y 

+ {dtP^) log a 

+ 2 an^F*' [ 7 * - £>* logo - 

-2aKhlF’^ - {an.F’') - {dt - P^du) (KF^) 


( 2 . 8 ) 




S=-c 


K.jK'-^ - AKn^F'' - 2 [n^F'^ f - 2K^ 


(2.9) 


+ [dt - P'^du) {n.F’') + B, 


where 'yfj are the Christoffel symbols associated with Di, 7 ^ = Kij = {l/2a){dt — Fp)hij is the extrinsic 

curvature of the Cauchy foliation and the quantities B,B^,B^j are constraint modification terms defined in (2.20)- 


(2.22) which are used to stabilize the constraint propagation system. With the addition of constraint modification, 
(2.4)-(2.6) correspond to equations (2.33), (2.36) and (2.38) of [TJ, as corrected for typographical errors. See [IH] for 
further details. 

Constraint preservation for the Cauchy problem follows from the system of homogeneous wave equations satisfied 
by the harmonic constraints (2.3). Constraint propagation can be extended to the initial-boundary value problem 


by implementing the hierarchy of Sommerfeld outer boundary conditions presented in I1I2D]. However, this has not 
yet been incorporated in the version of the code being used here and we maintain constraint preservation by causally 
isolating the region of interest from the outer boundary. 

For the purpose of using the method of lines to apply a Runge-Kutta time integrator, we re-write the system 


(|2.4|)-(|2.6|) as first differential order in time and second differential order in space by using the timelike vector field 

( 2 . 10 ) 


1 


n'' = - ((50 - wP'^) 
a 


to introduce the auxiliary variables 


n := = h^d^a 

n* := = fi^d^p^ - p’^dkF 

Bij := '^Fagij = — [fi^d^^hij -|- 


( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 


Here the function w is chosen to be unity everywhere except near the boundary, where it smoothly goes to zero. 
Insertion of (2.10) into (2.II)-(2.I3) gives the evolution equations for the lapse, shift and 3-metric, 


dta = aH -|- wP^dta 

dtp^ = an* + -p^P^dka - P^p'^dkW 

a 

dthtj — 2aTltj a dj^htj gitidjfi^ -t- g^jdtfi^'^ 


(2.14) 

(2.15) 

(2.16) 


The evolution equations for the auxiliary H-variables then follow from the first time derivatives of (2.14)-(2.16), after 


using (2.14)-(2.16) and (2.4)-(2.6) to eliminate first and second time-derivatives of the metric variables. 


We modify the evolution system by terms vanishing modulo the constraints based upon the results presented in [3] . 
For this purpose we set 


A'**' = - 


ai 

y=5 




S + Cpa-C^C^ 


\/^ 


(2.17) 
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where e is a small positive number, are positive parameters and 




is a Riemannien 4-metric. In terms of 


^ i29p.f9p(y-^^ 


the B-teims added to (2.71-(2.9) for constraint modification are 


B = - (-/ 

a 




Bij — —2Aij. 


(2.18) 

(2.19) 


( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 


B. Numerical implementation 


The evolution algorithm uses fourth order centered finite differencing and sixth order Kreiss-Oiliger type dissipation 
on a cubic Cartesian grid. The outer boundary points are updated using a Sommerfeld algorithm provided by the 
Cactus toolkit [S] , while we use summation-by-parts finite difference operators to update the neighboring three points 
along the normal Cartesian axis. 

In addition, following the lead of the BSSN method for treating punctures, we compute the spatial derivatives of 
the 3-metric hij using the conformal rescaling 


hij = (2.23) 

to obtain 

dkhij = hijdkh^^^ + (2.24) 

dkdthij = hijdkdth^^^ -|- (dkhij^ (2.25) 

-I- (dkh^^^^ (dehij'j + h^^^dkdjiij, 


where derivatives of hij are computed using finite difference operators. The derivatives of the conformal factor 
are computed using 


dkh^'^ = h^'^dk 

log(hl/3) 





dkd,h^/^ = h^/^dk 

log(hi/3) 

di 


+ h^/^dkdi 



(2.26) 

(2.27) 


where the finite difference operators are applied to the quantity [log (^^^^)] ■ 

In addition to the Cactus infratructure we use Carpet mesh refinement [2TJ [22] . The MOTSs are tracked using the 
horizon finder developed by J. Thornburg [23||24], which decomposes a topologically spherical surface into 6 abutting 
cubical patches. This finder is based upon a search algorithm which assumes a star-shaped domain inside the horizon. 

As this is the first example of harmonic evolution of black holes using the puncture method, as opposed to excision, 
we adopt a conservative approach to insure that the puncture never lies exactly on a grid point. This is made possible 
in the simulation of a head-on collision by offsetting the puncture a half grid-step from the axis along which the 
collision proceeds. In addition, for simulation of a head-on collision, harmonic gauge forcing is not necessary and we 
set = 0. Further studies and code development would be necessary to better understand the harmonic application 
of the puncture method so that the code could handle generic binary inspirals. With the present techniques, however, 
we can evolve the head-on collision for time scales well past the formation of a common apparent horizon and are able 
to track the subsequent penetration of the individual MOTSs. 













6 


III. THEORY OF THE UNIQUENESS AND STABILITY OF MOTSS 


Again consider a 3+1 foliation x°‘ = (t, x*) with metric decomposition 


(3.1) 


where is the unit future directed timelike normal to the foliation. In a slice Air '■= {t = t} consider a 2-surface S 
which is the boundary of a set Q. We define its outward unit spacelike normal Ni to point out of O. If S is defined 
as the level set of a function s, then 


Ni = -^d^s 

Vd 


(3.2) 


with 


(3.3) 

(3.4) 


D = h^\dks){dis). 

Setting = h^^Ni (so that = 0), this leads to the further decomposition 

+ s^ij. 

The outgoing null direction normal to S is 

=n^ + N^^, (3.5) 

where the normalization is determined by the Cauchy slicing. The expansion in this outgoing null direction is 


e+ = P + H (3.6) 

where 

H = (3.7) 

is the mean curvature of S in the Cauchy slicing, 

P = (3.8) 

is the 2-trace of the extrinsic curvature of the Cauchy slicing and V is the covariant derivative with respect to g. 

5 is a MOTS if 9+ = 0, i.e. P + H = 0. Similarly, 


e_=P-H (3.9) 

is the expansion in the ingoing null direction 

= (3.10) 


normal to S. The MOTS S is inner trapped if 6*_ < 0. Since 9- = 9+ — 2H = —2H for a MOTS, the trapping 
condition is equivalent to H >0. 


A. Stable and outermost MOTSs 

For a given slice Mr and a MOTS S C Mr one can consider the normal graphs of a function u G C°°(S), i.e. 
the surface parametrized by 


Fu : S —)■ Air : P '->■ exp(uN^) (3-11) 

where exp denotes the exponential map of Air- The operator 0+ : C°°(S) —>■ C°°(S) that maps a function u to the 
pull-back of the value of 0+ on to S has linearization given by 

Lf = -A/ - 2s^^SAdBf + fi^Rs - SaS^ + D^Sa - Ix^^Xab - 


(3.12) 
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with A the surface Laplacian of S, Rs denotes the scalar curvature of 5, Xtiv = denotes 

the Einstein tensor of the space-time and capital letters refer to intrinsic coordinates for S. In particular 
denotes the inverse of the tangential projection of to S. 

Although L is not self-adjoint, there exists a real eigenvalue A(iS) G cr(L), which is the unique minimizer for the 
real part in the spectrum a{L) of L. This eigenvalue is simple, and the corresponding eigenfunction </> can be chosen 
with a definite sign. If A(5) > 0 we say that S is stable, if A(5) > 0 it is called strictly stable. We refer to [HI [15] for 
further details. 

For the following, we require the surfaces S in question to be either bounding, S = dfl, or bounding with respect 
to an interior boundary dM., that is 5 = dft \ dM.. In both cases, we write S = S’*"!! or refer to S being an outer 
boundary in this situation. In the scenario considered here, the inner boundary exists and is formed by trapped 
surfaces enclosing the punctures. 

A MOTS S = S+n is called outermost if for all other MOTSs S' = 9+II' with II' D O it follows that H' = H. In 
other words, there are no MOTSs on the outside of S. 

In [25] it was shown that if M.t contains bounding outer trapped surfaces, as is the case if Mr is asymptotically flat, 
then there exists an outermost MOTS 5°“* that bounds the trapped region in Ad,-. This means, it is the enclosure of 
the region that contains outer trapped surfaces. Note that 5°“* is not necessarily connected. All components of 
are stable, and 5°“* has area bounded uniformly from above by a constant depending only on the geometry of the 
slice Mt- Furthermore, it has the property that there exists a positive i5 > 0, again depending only on the geometry 
of the slice Air, such that any geodesic starting on 5°^^* in direction of its outward normal A* does not intersect 
within distance <5. In particular, two distinct components of have distance at least <5. The constants mentioned 
depend in particular on an intrinsic curvature bound and on bounds for the second fundamental form V and its 
derivatives. For the details of these estimates and all the dependencies of the constants we refer to [5^ . Similar results 
have also been derived in ESI [17]. In the situation considered here, Galloway [H] established that 5°“* is a union of 
topological spheres. In [IS] it has been furthermore established that the second fundamental form ViAj of a stable 
MOTS S is bounded in the supremum norm, provided the geometry of the slice is bounded. 


B. The maximum principle for MOTSs 

A useful tool in the analysis of MOTSs is the strong maximum principle. To state it, assume that for a = 1, 2 
are two connected G^-surfaces with outer normals Assume further that there is a point p such that and ^2 
touch at p. If the outer normals agree, Nl = N^, at p and S 2 lies to the outside of iSi, that is in the direction of Nl, 
and furthermore 


supd+[iSi] < inf 0_i_[iS2] (3.13) 

5i S 2 

then = ^ 2 . This version can be found in |25[ Proposition 2.4] or in |16j . It implies in particular that two distinct 
MOTSs and ^2 can not touch in such a way that their normals point in the same direction and one is enclosing 
the other. 

The strong maximum principle provides an interpretation of strict stability. Assume that 5 is a strictly stable MOTS 
with outward normal A®. Let be the principal eigenfunction of L. Deforming S in direction of the vector field (/>A* 
then yields a foliation of a tubular neighborhood U oi S with the following properties. First U\S = U~ U U~^ where 
A® points into U~^. Moreover, is foliated by surfaces with 0+ > 0 and U~ is foliated by surfaces with 0+ < 0. The 
maximum principle then implies that there are no surfaces S' C which bound relative to S and have d+[iS'] < 0. 
Furthermore, there is no MOTS in U with outward normal aligned with that of S. 


C. Evolution of MOTSs to MOTTs 

Regarding the evolution of MOTSs there are different approaches. In m it was shown that a strictly stable 
MOTS S can locally be continued to a smooth space-time track of MOTSs, i.e. a marginally outer trapped tube. 
More precisely, for a given f such that A4f contains a strictly stable MOTS Sf there exists e > 0 such that for all 
T G (f — s,f + s) there is a stable MOTS Sr in Air such that the Sr form a smooth space-like manifold. To emphasize 
the role of the stability operator in this picture, we recall the argument from El- Assume that Sr is a smooth family 
of MOTSs passing through Sf ■ Then we can parametrize this tube by a map 

: (f — £, f + e) X Sf ^ [J Mf 

TG(T — €,r-i-s) 


(3.14) 


such that = V^, where is perpendicular to Sr at each point along the tube. Note that has the decompo¬ 
sition 


= an^ + = a{n^ + N^) -h (7 - a)N^ = -h fN^, (3.15) 

where as before a denotes the lapse function of the slicing. Calculating the change of 0+ under the deformation by 
at time t, we thus obtain 


— ^aii^Sj^\Sr\ -t 6f]\[i^9^[Sf] — —aW + Lf, 
where the first contribution is calculated via the Raychaudhuri equation with 

w = x^^xAB + G^jf^r, 


(3.16) 


(3.17) 


and the second part is just the definition of the stability operator (3.12). Since is tangent to a MOTT we have 
6vi^0+[Sr\ = 0 and thus 


Lf = aW. 


(3.18) 


The operator L and the function W are given by the geometry of Sr and the space-time geometry of Air, whereas 
/ is a function determined by (3.16). If Sr is strictly stable then L is invertible. Thus the previous calculation can 
be turned around to conclude the existence of the desired MOTT. The causal structure of the tube follows b y the 
observation that the null energy condition implies IT > 0. From stability one then obtains / > 0 via equation (3.18) 
as in [Ml Lemma 3]. 

A similar argument was used in m to construct a MOTT through Sf in the case it is stable but not strictly stable 
assuming IT > 0 and IT ^ 0. 

A different approach is to take the outermost MOTS 5°“* of all slices Air and define the apparent horizon of the 
slicing as 71 = In generic space-times H is smooth up to a discrete set of outward jumps. See [13] for the 

details and the particular notion of genericity used therein. The same reference also provides the causal character of 
71, namely that 71 is achronal. Moreover, if Or denotes the interior of 5™*, then J+(Or) H Alo- C Lla for all t < a. 
Here J“''(Or) denotes the causal future of Or. In particular, if there is a MOTS at an initial time f then it will persist 
for all times. 


D. Approaching MOTSs 

Assume that the space-time and the slicing are completely regular. Then the constant <5, the bound for the area, and 
the length of the second fundamental form of 5°^^* remains uniformly controlled. If in this situation two components 
of a MOTS are closer than 5, neither of these two components can be part of 5°“*. In particular, the evolution of 
^out discontinuous at some stage of the evolution. By the causality of 71, the jump is outward, and a new 

outermost MOTS has formed outside the tubes of the two original ones. Generically, after the jump time the jump 
target will split into two branches of MOTTs, a stable branch traveling outward and an unstable branch traveling 
inward. For a complete discussion of this jump in the outermost MOTS, we refer to m- An important fact to point 
out is that if two MOTSs are close to each other then there is no reason to expect that they be stable, in contrast to 
the stability of 

Note that it is not known whether the area of 5°“* is monotonic across the jumps. For MOTSs in general, including 
it is not even clear whether monotonicity of area should be expected along smooth pieces. However, if the 
MOTS has positive mean curvature, so that H = —P > 0, then its expansion in any outward spacelike direction 
+ an^, where 0 < a < 1, must be positive. As a consequence, the continuous portion of the MOTT traced out by 
a stable MOTS with H > 0 must have a monotonically increasing area. Such MOTSs will also be inner trapped, i.e. 
6 »_ = P - 77 < 0. 


E. Exterior osculation of MOTSs 


The following theoretical observations pertain to the collision between the MOTSs of a large and small black hole 
as they first touch. Note that this situation is not prevented by the strong maximum principle. However, at the time 
of first contact, a common horizon already has formed according to the previous section HI D| 
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In (3.6), the contribution to 0+ from P is common to both MOTSs, so that at a common point of osculation we 
must have 


IT _ Tpr 

(small) (large) ’ 


(3.19) 


Here the two mean curvatures are defined with respect to the respective outer normals, which in this case have 
opposite orientations. Thus the two MOTSs do not share the same outgoing null direction i. As a result of (3.19), 


the mean extrinsic curvature radius of the large and small MOTSs must match at the point of osculation. This can 
be deceiving in terms of a coordinate picture of the MOTSs since the connection of the Cauchy slice enters into the 
mean curvature. We have 


H = s^^d.Nj - s^^VtNk. 


(3.20) 


At the point of osculation Afi(sma«) = ~ ^i(ia rae) so that the second term has the same magnitude but opposite sign 
for the small and large MOTSs. Hence (3.19) implies 




(3.21) 


or 


^(small)^ ^ ^ (large)^ ^ ^ (<■<>■‘^9^) N^^iarge) ■ (3.22) 

This relates the coordinate curvatures didjs of the functions describing the MOTSs, as provided graphically by the 
code output. More geometrically meaningful output are plots of H during the evolution of the large and small 
MOTSs and plots of the time dependence of their surface area, at a sequence of times elapsed during the evolution. 
An important property is whether H > 0 so that the MOTSs are trapped. 

In Thornburgh’s apparent horizon finder [231124) . s = r— h{y^), where r is a standard radial coordinate measuring 
Euclidean distance from some point Xq and are spherical coordinates arising from a six-patch treatment of the 
unit sphere. Then 


diS = 


X *' _ ryt 

U/f) 


- dih. 


(3.23) 


In the axisymmetric case corresponding to the head-on collision of black holes, we must have dih = 0 on the symmetry 
axis and therefore also at the points where an osculation can occur. Thus, at an an osculation point, we have 
diS = (cc® — XQ)r~^ and D = hij{x'^ — Xq)(x-^ — Xn)r~^ = hxx, where we align the symmetry about the cc-axis. Thus 


D(^smaii) = D(^iarge) ^t the point of osculation and (3.22) reduces to 


s^^didjS(^ 


ill) ^ ^(large) 2P^r>l^dkS(^large), 


(3.24) 


where d^si^iaj-g^) — ^fc'5(sma/q- 


F. Interior Osculation 

If the MOTS of the small black hole were to remain intact and continue to completely penetrate the large one 
then there would again be a point of osculation between them just as the small one completely passes inside. At this 
contact point p, the two outer null directions have the same orientation. Thus 0+ (p) = 0 for the common outgoing null 
direction i{p). The strong maximum principle implies in this case that the two MOTSs coincide globally. Thus the 
mean curvatures of the large and small MOTSs must adjust dramatically to match each other unless full penetration 
is obstructed by a singularity or some other feature. One such obstruction of an artificial nature can arise due to 
singularity excision. Here we treat the singularities by punctures, which can restrict either the simulation or the 
horizon tracker to the point of half penetration. 

According to section [IH B[ it is impossible in this case that either of the two osculating MOTSs is strictly stable 
up to the time of coincidence. In fact, the normal separation of one MOTS above the other yields after linearization 
at the time where they coincide a function in the kernel of the stability operator. If the approach is fast enough, this 
function is non-zero and changes sign. This implies instability of the two individual MOTSs shortly before and at the 
time of coincidence. 

If the two individual inner MOTSs do coalesce into a single one, this raises the possibility of a scenario in which 
the two MOTTs traced out by the individual MOTSs merge and the resulting MOTT connects to the unstable 
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FIG. 1. A possible scenario for the evolution of the trapped tubes. The outermost MOTSs are drawn in bold. The common 
outermost MOTS develops before the first time of contact. After the separate MOTTs penetrate, the diagram illustrates the 
speculative scenario that they merge and join with the unstable branch from the jump. This speculative part of the figure is 
indicated by the dashed part . 


branch of the common horizon, whose stable branch is the MOTT traced out by the apparent horizon, as described 
in section [III D and depicted in Fig. The figure illustrates how the two individual MOTTs can coalesce and join 
continuously to the outer apparent horizon. Although this figure looks highly suggestive, there is no reason to believe 
that it is actually valid. What currently is known is that there always exists an outermost MOTS, that it can jump, 
and that at the time of the jump there are two branches emanating from the jump target. The current state of the 
theory is not able to determine how long the unstable MOTSs continue to exist. See m for a related discussion. 


IV. SIMULATIONS 

Here we present the results of the simulation of the head-on collision of two black holes with initial masses m and 
m/4. Besides outputting the coordinate shapes of the MOTSs detected by the horizon finder, in order to analyze 
their intrinsic geometrical structure we also output plots of the time dependence of their mean curvature H and total 
area. Since this is the first application of the harmonic code described in Sec. |llj we also present convergence tests 
to establish validity. Convergence of the constraints is checked in the region of interest close to the MOTSs. This 
avoids the expected loss of convergence near the punctures and near the outer boundary (where constraint preserving 
boundary conditions have not been implemented). Because the run times necessary for the head-on collision are short, 
errors from the outer boundary treatment do not effect the region of interest. 


A. Initial configuration of the MOTSs 

We use the same time symmetric Brill-Lindquist data to initiate an axisymmetric head-on-collision as in with 
the exception that the black holes now have unequal masses. The conformal flatness of the initial 3-metric provides 
a natural choice of Euclidean coordinates {x,y,z). The initial “bare” masses of the punctures, neglecting their 
interaction, are mi^iarge) = ‘^fn^smaii) = 0.8M, corresponding to the Euclidean coordinate radii of the non-interacting 
black holes = QAM. We choose the a;-axis to be the axis of symmetry, with the punctures 

corresponding to the two black holes initially separated by Aa; = l.OM. Due to their interaction at this separation, 
the actual masses of the individual black have ratio M (large) ~ (small)- (See [l2] for a discussion of the leading 

order effect of a finite separation on the physical parameters.) In order to avoid numerical problems, the punctures 
are offset by a half grid-step from the axis of symmetry. 

This initial configuration reduces the initial distortion of the black holes while also reducing the time required 
for merger. The sequence of three time symmetric initial data sets in Fig. illustrates the distortion that would 
arise from smaller separations. It also illustrates some of the problems which arise in interpreting the output of the 
numerical simulations and emphasizes the importance of monitoring geometric invariants in tracking the evolution of 
the MOTSs. 
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FIG. 2. Cross-sections of the coordinate shapes, in units of M, of the individual MOTSs which approach each other in a 
sequence of time-symmetric Cauchy slices. The uniqueness property of minimal surfaces keeps the MOTSs from touching and 
creates considerable coordinate distortion when they are close. 


At a moment of time symmetry the MOTS condition reduces to i? = 0, i.e. from (3.20) 


Then for a sequence of initial data for which the two black holes approach each other, (3.24) implies 

S ^ y ^ ^ (large) 


(4.1) 

(4.2) 


at the point of closest approach. Thus if the coordinate shape of one MOTS appears convex at the point of closest 
approach then the other must appear concave with the exact same magnitude of curvature. We see this effect in the 
initial data sets shown in Fig[^ 

Another feature of the sequence of time symmetric data sets is that the MOTSs cannot touch as their initial 
separation is made smaller. This is also a consequence of the uniqueness theorem for MOTSs. In the time symmetric 
case both MOTSs satisfy 9+ = 0- — 0. If they were to touch then at their common point the outer null normal to 
one MOTS would be the same as the inner null normal to the other. But since both null directions have vanishing 
expansion, this would violate the uniqueness theorem. Note that in the time symmetric case, a MOTS is also a 
minimal surface so that this result also follows from the uniqueness theorem for minimal surfaces. 


B. Approaching MOTSs 

As the individual MOTSs approach a common outer horizon (apparent horizon) forms at f = 0.144M, which is 
before they touch in accord with the theory described in Sec. |III D[ Figure shows the individual MOTSs at time 
t = 0.384M and at t = 1.452M just after the appearance of the common horizon. At t = 1.452M, the bifurcation of 
the common horizon has produced a stable branch propagating outward and an unstable branch propagating inward. 
Up to this stage, except for the asymmetry produced by the unequal black hole masses, the results are qualitatively 
similar to those found in El- 

Surface plots of the mean curvature H of the individual MOTSs, on the 6-patch spherical coordinates of the horizon 
finder, at times t = 0.384M and at t = 1.452A/ are shown in Figure]^ Initially, these mean curvatures are zero, as a 
result of the time symmetry. At t = 0.384M, they are fairly uniform, with H/^sraaii) ~ 2 and H(iarge) ~ 0-5 in roughly 
the 4 to 1 ratio of the initial masses. At t = 1.452M, the larger MOTS has undergone signihcant distortion due to 
the smaller one, although this is not evident in its coordinate shape. 


C. External osculation 


The front-sides of the approaching MOTSs touch att = 1.97384M, as portrayed by their coordinate representations 
in Fig.[^ At this time, Fig.|^shows surface plots of the mean curvatures of the individual MOTSs. Now the distortion 
has increased so that their mean curvatures are equal at their common point, as required by (3.19). 
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FIG. 3. Coordinate shapes of the individual MOTSs at times t = 0.384M (left) and t = 1.452M (right). At the earlier time, 
the coordinate shapes are almost spherical and then begin to distort as the MOTSs approach. 



FIG. 4. The mean curvature H of the individual MOTSs at times t = 0.384:M (left) and t = 1.452M (right) laid out 
on the angular coordinates p and a of the six cubed-sphere patches of the horizon finder. At the earlier time, the mean 
curvatures are fairly constant and have roughly adapted from their initial time-symmetric values iF = 0 to the expected values 
H(smaii) ~ '^H^iarge) for the Corresponding masses. On the right, considerable distortion has formed on the front end of the 
large MOTS. 


D. Penetrating MOTSs 

After osculation, the individual MOTSs continue to penetrate. Figure shows the coordinate shapes of both 
MOTSs as the penetration proceeds. The left panel of Fig. [^shows the two individual MOTSs at t = 2.1785M. The 
coordinate shape of the larger MOTS is clearly being affected along the axis where the small one has penetrated. 
The right panel of Fig. [^illustrates the coordinate shapes at a later time t = 2.538M, when the small MOTS has 
penetrated approximately half-way inside the larger one. At this stage, the puncture near the center of the small 
MOTS can be expected to have a signihcant effect on the larger one. This is evident in Fig.j^from the sharp localized 
feature that has developed at the front of the large MOTS. 

As the simulation progresses further, the horizon finder is not able to locate the larger individual MOTS. This 
could be due either to the horizon hnder breaking down for numerical reasons due to the puncture or simply due to 
the larger MOTS ceasing to exist. Further code development would be necessary to resolve this issue. 

The left and right panels of Fig. give surface plots of the mean-curvature of the individual MOTSs at the same 
times t = 2.1785M and t = 2.538M as in Fig. The left panel of Fig. shows that the mean curvature at the front 
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FIG. 5. Coordinate shapes of the individual MOTSs at time t = 1.97384M when they touch. The picture on the right zooms 
in on the contact point, at which the two MOTS have equal mean curvature. 



FIG. 6. Plots of the mean curvatures H of the individual MOTSs laid out on the (p, a) coordinates of the six patches of the 
horizon finder at the time of osculation t = 1.97384M. Their mean curvatures now agree at the contact point, as required by 
the theory. This is mainly achieved through the distortion of the large MOTS. 


of the larger MOTS has continued to grow and now overshoots the mean curvature of the smaller MOTS. At the later 
time displayed in the right panel of Fig. the mean curvature at the front of the larger MOTS has now fallen back 
below the mean curvature of the smaller MOTS and has developed some highly oscillatory behavior. Again, this could 
be the result of the puncture inside the smaller MOTS as it crosses the surface of the larger MOTS at t « 2.45M. 

Figure|^shows the time evolution of the mean curvatures at the (facing) front ends of the individual MOTSs as they 
approach and penetrate. These are the points around which their mutual distortion is greatest. The mean curvature 
at the front end of the smaller MOTS grows steadily from its initial time symmetric value H = 0 and then settles to 
a fairly constant value as its geometry adjusts to the ambient gauge. In contrast, the mean curvature F7(iai.ge) sit the 
front end of the larger MOTS shows a rapid increase as the smaller MOTS approaches and penetrates. At the time of 
osculation, t = 1.97384M, both mean curvatures agree, as required by (3.191. As the penetration continues, H(j^arge) 
continues to rise until the puncture inside the small MOTS crosses it. 7I(/arge) then first drops discontinuously to a 
much smaller value and then begins to rise again during the period when both MOTS overlap. Runs with different 
spatial resolutions and different time steps confirm that this jump is not a numerical artifact of the evolution code. It 
could possibly result from a discontinuity in the MOTT trace out by the large MOTS similar to the jump that occurs 
in the formation of the outer common horizon. Or it could be an indication of a more complicated horizon structure 
which is difficult for the horizon finder to fully describe. This is another issue that we defer to future work. 
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FIG. 7. Zoomed-in view of the coordinate shapes of both individual MOTSs during penetration at times t = 2.1785M (left) 
and t = 2.538M (right). The puncture near the center of the small MOTS has produced a sharp effect on the front end of the 
large MOTS. 



FIG. 8. Mean curvatures H of the individual MOTSs laid out on the (p, a) coordinates of the of the horizon finder at times 
t = 2.1785M (left) and t = 2.538M (right). On the left, as the penetration first proceeds, the mean curvature at the front of 
the larger MOTS rises above that of the smaller MOTS. Later, on the right, the puncture of the smaller MOTS has had an 
oscillatory effect on the the mean curvature of the larger one. Gomparing the two panels, the mean curvature of the smaller 
MOTS has undergone little change. Further evolution becomes problematic beyond this stage. 


This discontinuous drop in H(jarge) has a counterpart in Fig. 10 which plots the time evolution of the areas A of 


both individual MOTSs. The larger MOTS shows a discontinuous jump in area at the same time its mean curvature 
undergoes a discontinuous drop. The small MOTS shows very small growth in area during the entire run, with an 
approximate range 3.941 < A(smaH) < 3.948. (Near the outset, A(^srnaU) undergoes a small anomalous decrease which 
can be attributed to numerical error since this would otherwise be inconsistent with the strict stability of the small 
MOTS at this time.) At the end of the run, the ratio A(^ii^j.ge)/A(^s 7 naii) has not undergone substantial change, which 
offers no clue regarding their possible future coalescence. 


E. Convergence 


In order to verify the new geometrical features we have found, we have measured the convergence of the harmonic 
constraints i.e. (2.31 with the gauging forcing = 0 used in our simulations. We monitor the constraints via 
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FIG. 9. Mean curvature H at the front ends of both MOTSs as a function of time. grows steadily from its initial time 

symmetric value and then settles to a fairly constant value. H^arge) rises rapidly as the MOTSs osculate and penetrate and 
then drops discontinuously after the puncture inside the small MOTS passes through. 
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FIG. 10. Areas A of the small and large MOTSs as a function of time. Afiarge) undergoes a discontinuous increase at the same 
time that its mean curvature drops discontinuously. Afterwards, A^iarge) drops continuously but at a rapid rate. A^smaii), 
whose value is multiplied by 10 to match the scale of the plot, shows very little variation during the entire run, which offers no 
clue regarding the possible future coalescence of the two MOTSs. 


the norm 


|Cp := 


(4.3) 


where is the positive-definite metric defined in ( |2.18 ). Although the code has fourth order finite difference accuracy, 
the time integrator used in conjunction with the adaptive mesh refinement limits the overall accuracy to second order. 
Confirmation of second order convergence of \C\ in a region including the penetration of the MOTSs supplies the 
necessary code verification. 


Figure[^plots \C\ along the cc-axis, in the relevant region where the puncture and outer boundary are not included, 
for three different gridsizes, h = (3.75 • 10“^M,4.375 • 10“^M,5 • 10“^M), at time t = 2.45M. At this time, the 
individual MOTSs have penetrated and overlap in the region between x = 0.45M and x = 0.63M, as indicated by 
the vertical lines in the figure marking the front and back end of both MOTSs. The curves are rescaled in accordance 
with the expected second order convergence of the constraints. Second-order convergence is confirmed even very close 
to the punctures in the region where the MOTSs overlap. This confirms the validity of the geometrical features we 
have observed of the MOTSs and their penetration. 
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FIG. 11. Rescaled norms of the harmonic constraints \C\ along the x-axis for 3 different resolutions (r0,rl,r2), corresponding 
to the gridsizes h = (3.75 • 4.375 • 10“^M, 5 • 10“^M), at time t = 2.45M. The norms have been rescaled for resolution 

to confirm that the convergence is second order. The vertical lines bound the region that the penetrating MOTS overlap. 


V. CONCLUSION 


The results presented here give independent confirmation of the penetration of individual MOTSs reported in [T] 
and, by using punctures rather than excision, we are able to follow the penetration to approximately the halfway stage. 
While the current code also works well with the excision module developed for the PITT Abigel code (and later adapted 
to the AEI Harmonic code [T]), the excised domain restricts the simulation to a very small penetration. The horizon 
finder algorithm would fail as soon as its required search domain intersected an excised region. Unfortunately, even 
with the use of punctures, the tracking of the larger MOTS breaks down before the most interesting final stage. Either 
the individual MOTSs eventually cease to exist or hit the final singularity; or otherwise they must coalesce. 

The apparently discontinuous behavior of the large MOTS near the end of the simulation indicates that there might 
be underlying complexities introduced by the time slicing. It is particularly interesting that the jump leads first to 
an increase in its area, followed by a continuous fall. In order to investigate this feature and track the penetration 
further, it might be necessary to redesign the horizon finder or to make some fortuitous choice of harmonic gauge 
source function. This is an open problem since this is the first binary black hole simulation carried out in the harmonic 
formulation using punctures. 

One attractive approach to the code modifications necessary to continue our simulations is to fill the puncture 
with artihcial matter, using a version of the “turduckening” procedure [30) . This is feasible because the uniqueness 
theorem governing the coalescence is independent of a positive energy condition on the matter. 

The penetration of MOTSs is a new aspect of the of the dynamics of MOTTs. The possibility exists that the 
individual MOTTs merge and connect back to the (unstable) branch of the common outer horizon to form a globally 
connected MOTT, This makes the system of penetrating MOTSs a rich laboratory for the geometry governed by 
Einstein’s equations. 
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